Discontinuous nonequilibrium phase transitions in a nonlinearly pulse-coupled 

excitable lattice model 
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We study a modified version of the stochastic susceptible-infected-refractory-susceptible (SIRS) 
model by employing a nonlinear (exponential) reinforcement in the contagion rate and no diffusion. 
We run simulations for complete and random graphs as well as d-dimensional hypercubic lattices 
(for d = 3,2, f). For weak nonlinearity, a continuous nonequilibrium phase transition between an 
absorbing and an active phase is obtained, such as in the usual stochastic SIRS model [Joo and 
Lebowitz, Phys. Rev. E 70, 036II4 (2004)]. However, for strong nonlinearity, the nonequilibrium 
transition between the two phases can be discontinuous for d > 2, which is confirmed by well- 
characterized hysteresis cycles and bistability. Analytical mean-field results correctly predict the 
overall structure of the phase diagram. Furthermore, contrary to what was observed in a model of 
phase-coupled stochastic oscillators with a similar nonlinearity in the coupling [Wood et at, Phys. 
Rev. Lett. 96, 145701 (2006)], we did not find a transition to a stable (partially) synchronized state 
in our nonlinearly pulse-coupled excitable elements. For long enough refractory times and high 
enough nonlinearity, however, the system can exhibit collective excitability and unstable stochastic 
oscillations. 

PACS numbers: 05.70.Ln,05.50.+q,02.50.Ga,87.15.Zg 



I. INTRODUCTION 

Understanding collective effects of noisy excitable ele- 
ments is essential for several disciplines, such as ncuro- 
scicncc, epidemiology, and chemistry, among others. An 
isolated excitable element is a dynamical system vifhich 
stays in a quiescent state until it suffers a sufficiently 
strong perturbation. In that case its trajectory in phase 
space can be characterized by an excited state, which 
is then followed by a refractoriness to further pertur- 
bations before returning to rest. A minimal (discrete) 
model of an excitable element therefore consists of a 
three-state system In the parlance of neuroscience 
(epidemics), each element could represent a neuron or 
patch of active membrane (individual) which sequentially 
becomes polarized (susceptible), spiking or depolarized 
(infected), and then refractory (recovered). Collective 
dynamics emerge because quiescent elements are typi- 
cally perturbed by excited elements. 

A simple system which incorporates these ingredients 
in a scenario with noise is the stochastic susceptible- 
infected-recovered-susceptible (SIRS) epidemic lattice 
model 0. It is a continuous-time model in which a site 
goes from susceptible to infected at a rate which depends 
on the density of its infected neighbors. In epidemiol- 
ogy, the model (or variants thereof) can be employed to 
investigate, e.g., whether an initial density of infected 
sites (which is usually chosen as the order parameter) 
will reach a nonzero stationary value or decrease to zero. 
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If all sites become quiescent, the dynamic halts and the 
system is said to be in a absorbing state Q . 

Increasing the coupling between infected and sus- 
ceptible sites, the SIRS model undergoes a continuous 
nonequilibrium phase transition from an absorbing to an 
active phase characterized by a stationary nonzero den- 
sity of infected sites. However, experimental data (from 
neuroscience, epidemiology, and chemistry, among oth- 
ers) can exhibit also global oscillations. This additional 
transition has been observed mostly in cellular automata, 
where the sites arc synchronously updated [l|, 0, H @1- 
This technical detail is apparently relevant since global 
oscillations in stochastic continuous- time models are less 
common. As discussed by Risau-Gusman and Abram- 
son, they sometimes appear as stochastic oscillations in 
single runs of the model but disappear in trajectories of 
the averaged lattice activity and analytical descriptions 
(usually mean- field) thereof 0. Global oscillations pre- 
dicted by mean-field (MF) approximations were observed 
in both non-Markovian Q and Markovian [§] models of 
three-state continuous-time stochastic oscillators. These 
models, however, do not have an absorbing state. 

Here we investigate a modified version of the SIRS 
model with a nonlinear rate which can be considered 
a pulse-coupled excitable version of the original phase- 
coupled oscillator model of Wood et al. @. We will show, 
on the one hand, that this nonlinear (but Markovian) ex- 
tension of the SIRS model is apparently insufficient for 
the generation of sustained collective oscillations. On the 
other hand, the model presents phase transitions into 
an absorbing state which can be continuous (if weakly 
nonlinear) or discontinuous (if nonlinear enough and for 
d > 2). Wc therefore provide a three-state continuous- 
time model which can undergo discontinuous phase tran- 
sitions like those of Bidaux et al. for a two-state cellular 
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FIG. 1; (Color online) Phase space [P{S), P{I)] showing the nuUclines P{I) = (dashed lines) and P{S) — (dot-dashed 
lines). Grey triangles corresponds to the physically acceptable region < P{I) < 1 — P{S)- In each row, a is kept constant and 
b increases from left to right. In the lower panel (b) (a = 0.5), the absorbing state at P{S) — 1 loses stability in a transcritical 
bifurcation as b increases. In the upper panel (a) (a = 2.25), the P{I) = nuUcline is bent by the exponential nonlinearity [see 
eq. 13] , creating an active state with finite P{I) in a saddle-node bifurcation. Closed (open) symbols denote stable (unstable) 
fixed points. Note that the nuUcline P{I) = always comprises the line P{I) — 0. When the absorbing state is unstable (a 
saddle), P{I) = corresponds to its stable manifold (rightmost column). Except for this case, the trajectories (solid lines) 
represent the stable and unstable manifolds of the saddles. 



automaton [10[ except that the nature of the transition 
here can be controlled by a free parameter (not only the 
spatial dimension [l^l)- Moreover, differently from most 
models presentin g n onequilibrium discontinuous phase 
transitions [1, [Til Il2| , the parameter which controls the 
nature of the transition is not diffusion (which is absent 
from our model). The model can also exhibit a differ- 
ent bifurcation scenario from what is usually observed 
in nonequilibrium lattice models, showing collective ex- 
citability and unstable global oscillations. 

The paper is organized as follows. In Sec.[TI]we intro- 
duce the model, analyze its mean-field solution, and com- 
pare it to simulations of random and complete graphs. 
Simulation results for d-dimensional hypercubic lattices 



arc shown in Sec. IIIIl while our concluding remarks are 
discussed in Sec. HV] 



II. MODEL 

In the conventional stochastic version of the SIRS epi- 
demic model, a susceptible (S) site at position x (x = 
1,...,N) becomes infected (/) at a rate Xnj{x)/k{x), 
where k{x) is the number of neighbors of x, out of which 
nj{x) are infected, and A is the so-called infection rate 
(the coupling parameter). After that, the site becomes 
temporarily insensitive to its surroundings (hence the 
term pulse coupling), jumping from infected to recov- 
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ered (R) at a constant rate S, and from recovered back 
to susceptible at a rate 7. This is summarized as follows: 

S — >I atrateAn//fc, (1) 

/ — > R at rate (5, (2) 

R — > S at rate 7. (3) 

Since all rates could in principle be normalized by 5, in 
the following we set ^ = 1 without loss of generality. 

For low values of A, an initial density of infected sites 
eventually dies out and the system reaches its unique ab- 
sorbing state (all sites susceptible). For A larger than a 
critical value Ac, on the other hand, a phase with nonzero 
density of infected sites becomes stable in the thermody- 
namic limit N ^ 00 (though the absorbing state is ob- 
viously always a solution). The transition at Ac (studied 
in detail by Joo and Lebowitz 0) is widely believed to 
be continuous and belonging to the directed percolation 
(DP) universality class [ij, . 

Adapting the nonlinear coupling employed by Wood et 
al. to pulse-coupled excitable elements, one obtains a 
generahzation of the SIRS model. Instead of Eq. [U we 
propose 

S — > / at rate 
g {ni/k, ns/k) = b [e<^i~^s)/k g-a«s/fc] ^ (4) 

where a and b are coupling parameters, is the num- 
ber of susceptible neighbors, and Eqs. [5] and [3] remain 
unaltered. Note that the interaction occurs only among 
first neighbors. The second term in rate |4] guarantees the 
existence of an absorbing state: if all sites are susceptible 
(thus nj = for all sites), they will remain susceptible 
forever. For small values of a (a <C 1), one recovers the 
linear behavior of the original SIRS model, with A ~ ab, 
to first order. So for small a, increasing b leads to a con- 
tinuous phase transition just like in the SIRS model 0. 
For large enough a, however, we will show that increasing 
b leads to a discontinuous phase transition. 

As a motivation to rate IH let us note that the dy- 
namics underlying neuronal firing is highly nonlinear in 
several of its aspects: membrane depolarizes when (typ- 
ically Na+) channels open very quickly, in a threshold- 
like behavior. This process, on its turn, is triggered by 
a (nonlinear) sum of smaller depolarizations induced at 
synapses by presynaptic neurons. Synaptic dynamics (in- 
cluding neurotransmitter binding, in the simplest case) 
is itself nonlinear. It is therefore not unusual that re- 
duced models of collective neuronal phenomena allow for 
nonlinear terms p^ . In our case, the nonlinearity is con- 
trolled by parameter a. 



A. Mean-field analysis 

The structure of the phase diagram in the (a, b) plane 
can be captured most easily by analyzing the mean-field 
version of the model. Letting P{a) be the probability 
that a site is in state a {a £ {S, I, R}), one obtains 

PiS) = -g[PiI),PiS)]PiS)+^PiR), (5) 
P{I) = g[P{I),P{S)]P{S)-SPiI), (6) 

P{R) = SP{I)-jP{R). (7) 

These equations are exact for a complete graph when 
fc = iV — > 00 and rid 00 with Ua/k = P{oi). In what 
follows, we will employ the stationary value of the density 
of infected sites P{I)* as the order parameter. 

By employing the normalization condition P{S) + 
P{I) + P{R) = 1, one can eliminate P{R) and study the 
resulting two-dimensional flow of the mean-field dynam- 
ics in the [P{S), P{I)] plane, as shown in Fig.[TJ For low 
values of a (lower panel of Fig. [T]), the absorbing state 
P{S) = 1 loses stability in a transcritical bifurcation, 
giving rise to a stable active state with a density of in- 
fected sites which increases continuously from P{I) = 0. 
For large enough a, the discontinuous character of the 
phase transition reveals itself in the mean-field equations 
through a saddle-node bifurcation (upper panel of Fig.[T|): 
an active state appears with nonzero P{I), while the ab- 
sorbing state remains stable. As usual, this bistability 
is regulated by the stable manifold of the saddle, which 
separates the basins of attraction of the two stable fixed 
points. Increasing b further, another transcritical bifur- 
cation occurs: the absorbing state loses stability and the 
active state becomes the only attractor of the system. 



B. Random graph simulations 

Hysteresis is one of the simplest fingerprints of mul- 
tistability, and in this system it can be clearly detected 
in simulations of Erdos-Renyi random graphs (l6j with 
finite average connectivity K, with which the mean-field 
equations show a good agreement. Simulations were per- 
formed for a fixed value of a. For each value of &, we al- 
lowed the system to evolve during tmax time steps. Each 
time step At (corresponding to N random updates Q) 
was chosen to be (i5 -I- 7 -f 6e°)~^ to make sure proba- 
bilities are less than one. Parameter b was increased or 
decreased in constant intervals 5b, and the initial condi- 
tion of the network for each value corresponded to the 
final condition of the preceding case. We applied a small 
rate h {h < l/N) to spontaneously excite quiescent sites, 
thus preventing the system from getting trapped in the 
absorbing state by finite-size fiuctuations [10, ■ 

The insets of Fig. [5] show that a loop in b leads to a 
hysteresis cycle [as observed for the density of active sites 
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FIG. 2: (Color online) Phase diagram of the (mean-field ver- 
sion of the) model. The dashed line marks the onset of phase 
bistability, with the sudden appearance of a stable active 
phase with finite order parameter [saddle-node bifurcation in 
Eqs.[5][7]- Solid lines denote a transition in which the absorb- 
ing state loses stability [transcritical bifurcation in Eqs. [5][7] . 
Symbols correspond to simulations in a random graph with 
7 = 1, if = 10, iV = 10", t^ax = 5 X 10^ and h = 10"^ aver- 
aged over five runs. The lower (upper) inset shows the order 
parameter P{I)* for a cyclic change in the coupling parameter 
b, showing a continuous (discontinuous) transition for a small 
(large) value of a. Upward (downward) triangles: increasing 
(decreasing) b. 
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P(/)] for high enough values of a (a > Oc). For low val- 
ues of a, the phase transition is continuous. The bound- 
aries between the phases in the (a, b) plane (main plot of 
Fig. [2) can be numerically obtained from the mean-field 
equations and easily estimated from the random graph 
simulations (hysteresis cycles did not change significantly 
when tmax was doubled; see also Sec. IIIip. 

To obtain the phase diagram shown in Fig. ^ we em- 
ployed a mean refractory time comparable to the excita- 
tion time, 7 = 1. However, we did not find qualitative 
differences in the mean-field phase diagram in the lim- 
iting case 7 — *■ oo, which suggests that the simulation 
results we present here will also be valid for a two-state 
system. 



C. Absence of sustained global oscillations 

In the process of scanning parameter space in search 
of collective oscillations, we have found small regions in 
which a Hopf bifurcation involving the active state can 
indeed occur. As we shall see, however, this does not 
necessarily imply the existence of sustained collective os- 
cillations, which we have not found in this model. 

We exemplify with results for 7 = 0.01, whose phase 
diagram is shown in Fig. [3^. Though qualitatively similar 
to Fig. [21 there is now a narrow interval of b values (for 



FIG. 3: (Color online) (a) Phase diagram of the (mean-field 
version of the) model for 7 = 0.01 (panels b-g also show 
single-run simulations of a complete graph with N = 10®). 
Bifurcations are as in Fig. [J] except for those additionally 
occurring within the marked rectangle (inset): for increas- 
ing b, a saddle-node bifurcation (dashed line) is followed by a 
Hopf (H) bifurcation (dotted line), after which a homoclinic 
(HC) bifurcation occurs (see text). The absorbing phase cor- 
responds to the grey region. The labels (b)-(d) in the inset 
show the points in parameter space which correspond to the 
phase portraits below, (b) A saddle and an unstable fixed 
point are born in a saddle-node bifurcation, (c) The fixed 
point becomes stable via a Hopf bifurcation and is surrounded 
by an unstable limit cycle, (d) After a homoclinic bifurcation, 
the limit cycle disappears. UUM = upper unstable manifold 
of the saddle; USM — upper stable manifold of the saddle; 
ULC = unstable limit cycle; LSM = lower stable manifold. 
The scale in (b) also applies (apart from a translation) to (c) 
and (d). The lower unstable manifold of the saddle always 
goes to the absorbing fixed point. Panels e, f, and g show 
time series for simulated trajectories, respectively, shown in 
panels b, c, and d, with examples of collective excitability and 
stochastic oscillations (see text for details). 



high values of a) in which the route to bistability, instead 
of the simple saddle- node scenario described in Fig.[lja), 
requires additional intermediate bifurcations psj [inset of 
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Fig. [11^ a)]. These occur while the absorbing fixed point 
remain stable. 

Starting from the absorbing phase and increasing 6, 
first a saddle-node bifurcation occurs in which the node 
is unstable [and quickly becomes a spiral, Fig. [Hb)]. At 
this stage, the system still has only a single stable fixed 
point, but the structure of phase space is such that the 
system has become collectively excitable: if the system 
is below the upper stable manifold (USM) of the saddle 
[Fig- ISHb)], it will monotonously return to the absorbing 
state, whereas a point above the USM will go through 
a long excursion around the upper unstable manifold 
(UUM) before coming to rest [Fig. [31[b)], displaying a 
spike- like time series [Fig. |3lje)]. 

Increasing b further, the spiral becomes stable in a Hopf 
bifurcation and is now surrounded by an unstable limit 
cycle (ULC) [Fig. [3Jc)]. Formally, in this parameter re- 
gion the system is collectively bistable, but note that the 
active phase will only be reached from initial conditions 
within the ULC [which, owing to the small value of 7, 
is also very small — see scale in Fig. EJb)]. Note also 
that the overcrowding of lines outside the ULC signals 
that it is weakly repulsive. The inner stable fixed point 
(the active phase) is correspondingly weakly attractive 
(which is the reason why the flux inside the ULC is not 
shown). This means that a collective oscillation solution 
exists (the ULC) but is so weakly unstable that it might 
get confounded with sustained oscillations [even in the 
numerical integration of eqs. [5][7]- From single-run sim- 
ulations of a complete graph with N = 10^, we have 
obtained a time series with oscillations [Fig. [3Uf)] which 
only disappear when, owing to finite-size fluctuations, 
the system reaches the absorbing state We there- 
fore confirm the scenario predicted by Risau-Gusman and 
Abramson: since the eigenvalues of the stable fixed point 
have a nonzero imaginary component, inevitable fluctu- 
ations will generate stochastic oscillations 0|- 

Finally, the ULC disappears in a homoclinic (HC) bi- 
furcation, after which the active fixed point is separated 
from the absorbing fixed point only by the stable mani- 
folds of the saddle. Note that the lower stable manifold 
[LSM, see Fig. [Hd)] no longer comes from the unsta- 
ble fixed point nor from the ULC but rather joins the 
USM as t -00. The LSM wiU gradually unfold as 
6 increases until a phase portrait similar to that of the 
central plot of Fig. [1^ is reached (before the saddle col- 
lides - for yet larger values of 6 - with the absorbing 
fixed point in a transcritical bifurcation). Collective ex- 
citability and stochastic oscillations remain present in 
this regime [Figs.[3Jd) andOJg)]. 

Our phase diagram emerged essentially from local sta- 
bility analysis, so in principle it does not exclude a saddle- 
node bifurcation of cycles from occurring within the ac- 
tive or bistable regions. However, numerical integration 
of Eqs. [5][7] for a variety of initial conditions and combi- 
nation of parameters did not show any signs of it. 

Although the above analysis is based on the mean-field 
approximation, it is worth mentioning that improvements 



on the mean-field approximation do not necessarily help 
the prediction of collective oscillations. Rozhnova and 
Nunes recently observed that the equations obtained 
by Joo and Lebowitz Q for the two-site approximation 
lead to sustained oscillations in a small region for 7 <C 1. 
However, when simulating random graphs in the same 
parameter region, these oscillations get damped [l9| . 



III. SIMULATIONS IN HYPERCUBIC 
LATTICES 

Since simulations with small values of 7 are very diffi- 
cult to perform , we now focus on the simpler bifurca- 
tion scenario of Fig. [5] and discuss the results of simula- 
tions for 7=1. Identifying the nature of the transition 
in hypercubic lattices is not so simple as for random and 
complete graphs. As has been recently discussed in de- 
tail by Takeuchi [l^, even a system which undergoes a 
continuous phase transition into an absorbing state (such 
as those belonging to the directed percolation universal- 
ity class) may show a hysteresis cycle when the coupling 
parameter loops around its critical value. This is due 
to the divergence of the transient times at criticality in 
the thermodynamic limit. In simulations, this gets re- 
flected in the width of the hysteresis cycle scaling with 
the ramp rate (defined as the increment in b per unit 
time) as {l/tmaxY^'''^'^^^ , where (3 is the critical expo- 
nent governing the order parameter [3, . 

Consider, for instance, the hysteresis cycles shown in 
Fig. [ll^a) for d = 3 and two values of t^ax differing by 
an order of magnitude (while the increments in the val- 
ues of b have been kept the same [l3|)- Whether or not 
the transition is continuous will depend on whether the 
width A6 of the hysteresis cycle shrinks to zero in the 
limit tmax 00. We have operationally defined Ab as 
follows: for each run, we attribute the upper end of the 
hysteresis cycle to the b value at which the density of 
active sites (averaged over tmax time steps) P{I) is first 
above I/^/N. Similarly, the lower end of the hystere- 
sis cycle is defined as the point where P{I) falls below 
The width of the hysteresis cycle A& is then ob- 
tained by averaging the difference over the runs. We 
plotted A6 versus 1/tmax in Fig. [4l[b) for two values of 
a. For a = 4.5 [which corresponds to the hysteresis cy- 
cles shown in Fig. H^a)] , a linear extrapolation leads to 
a nonzero value of A5 as tmax 00, which is consis- 
tent with a discontinuous phase transition. For d ~ 3 
and a = 1.875, on the other hand Ab decreases to zero 
as a power law (1/tmaxY''^^ , which is consistent with 
Takeuchi's prediction [l3| for the DP universality class 
[/? ~ 0.805(10) HI]. 

One could in principle feel uncomfortable with the 
above described criterion for deciding on the discontinu- 
ity of the transition because in practice the linear extrap- 
olation may not coincide with hmj^^^^^oo Ab if transient 
times are too long (and we expect them to be long near 
a continuous transition!). This problem becomes more 
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FIG. 4: (Color online) Simulations for d = 3 with iV = 25^ 
and h = 2.5 x 10~^. (a) (a = 4.5, averaged over 20 runs) 
The length of the hysteresis cycle A6 decreases for increasing 
traax [Ab dcuoted by the horizontal arrow for tmax = 10""). 
Upward (downward) triangles denote increasing (decreasing) 
b. Vertical lines show values of b employed in panels (c) and 
(d). (b) Dependence of Afo on tmax is qualitatively different 
for a above (filled circles) and below (open circles) ac, re- 
spectively, showing a nonzero or zero asymptotic value in the 
limit tmax — > 00. The lower plot is well fitted by a power 
law with exponent 0.55. Note that the precision in Ab is lim- 
ited by the increment in b employed in the hysteresis cycle 
(Sb = 1.8 X 10~^ in the case a = 1.875). Panels (c) and (d) 
show the time evolution of P{I) (averaged over 20 runs) for 
a = 4.5 and different initial conditions [with NP{I) sites in 
state / at t = and the remaining in state 5*] , showing phase 
bistability for b = 3.25 but not for 6 = 7. 



FIG. 5; (Color online) Simulations for d = 2 with N = 150^ 
and h = 2.5 x 10"''. (a) (a — 7, averaged over 10 runs) The 
length of the hysteresis cycle A6 decreases for increasing tmax 
{Ab denoted by the horizontal arrow for tmax — 2 x 10^). 
Upward (downward) triangles denote increasing (decreasing) 
b. Vertical lines show values of 6 employed in panels (c) and 
(d). (b) Dependence of A6 on tmax is qualitatively different 
for a above (filled circles) and below (open circles) a^, respec- 
tively, showing a nonzero or zero asymptotic value in the limit 
tmax 00. The lower plot is well fitted by a power law with 
exponent 0.63. Note that the precision in Ab is limited by the 
increment in b employed in the hysteresis cycle {Sb — 10~^ in 
the case a — 4.5). Panels (c) and (d) show the time evolution 
of P{I) (averaged over 20 runs) for a = 7 and different initial 
conditions [with NP{I) sites in state / at t = and the re- 
maining in state S], showing phase bistability for 6 — 4.4 but 
not for b = 21.5. 



salient as we decrease the spatial dimension, as depicted 
for d = 2 and a = 7 in Fig.O Note the smoothness of the 
largest hysteresis cycle in Fig. [Sl[a)jWhich is very similar 
to the ones observed by Takeuchi [l7[ near a continuous 
transition. According to the extrapolated value of Ab 
for tmax ^ oo, however, this transition would be consid- 
ered discontinuous [see upper plot in Fig. El^b)] , whereas 
for weaker nonlincarity [a = 4.5, lower plot of Fig. Elb)] 



we obtain again a power law {1/tmax)^'^^ [13] compatible 
with DP in d = 2 [/? ~ 0.583(4) [^]. 

To be sure of the discontinuity of the transition, it 
is simpler (and computationally less expensive) to in- 
vestigate directly the alleged bistability: we fix a and 
b and check the dependence of the stationary state on 
the initial condition [2l|, This test is shown in 

Figs. He), Hd), EJc), and [D^d) for = 3 and d = 2, 
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FIG. 6: (Color online) Phase diagram of the model in one-, 
two- and three-dimensional hypercubic lattices (squares, cir- 
cles and triangles, respectively). Open symbols (with dashed 
lines to guide the eye) mark the onset of phase bistability, 
with the sudden appearance of a stable active phase with fi- 
nite order parameter. Filled symbols (with solid lines to guide 
the eye) denote a transition in which the absorbing state loses 
stability. The region in between dashed and solid lines cor- 
respond to the bistable regime. Lines without symbols show 
MF results for comparison. Results correspond to simulations 
with 7 = 1 and /i = 2.5 x 10^^. For d = 1, 2, and 3, linear 
system sizes were L = 32400, L = 150, and L — 25 (where 
— L'^), the maximum values of tmax were 5 x 10^, 2 x 10''', 
and 10^, with results averaged over 10, 10, and 20 runs, re- 
spectively. 



respectively. Figures \^c) and EJc) confirm bistability 
since for lower (higher) initial values of P{I) the system 
converges to the absorbing (active) state. Figures IH^d) 
andinijd) serve as a negative control, confirming (in a re- 
gion where only the active state is stable) that the con- 
vergence to the absorbing state in the previous cases arc 
not due to finite-size fluctuations. We note that all sam- 
ples converged to the same attractor (either absorbing or 
active) as their average. 

The extrapolation limt^^^_,oo was employed to 
draw the phase diagram of the model for two- and three- 
dimensional lattices. As shown in Fig. [51 the qualita- 
tive structure of the phase diagram is well reproduced 
by the mean-field predictions, though quantitative agree- 
ment worsens as dimensionality decreases as expected. 
Note that the bistable phase for d = 2 is much smaller 
than for d = 3. For d = 1, the large error bars in Fig. [6] 
for large a emerge due to extremely large transients. Wc 
have not observed clearly discontinuous transitions for 
d — 1 up to a = 7. This is in agreement with the results 
of Bidaux et al. [l^, [2^ , as well as with Hinrichscn's con- 
jecture that discontinuous transitions in d = 1 should 
only occur in the presence of diffusion [23 |. 

IV. CONCLUSIONS AND PERSPECTIVES 

In summary, we have proposed a Markovian 
continuous- time lattice model of nonlinearly pulse- 



coupled excitable elements. Coupling depends linearly 
on rate b and nonlinearly on the dimensionless parame- 
ter a. We have shown that increasing the nonlinearity 
of the coupling leads to a change in the nature of the 
phase transition into an active state. In the regime of 
linear coupling (a <C 1), where the model approaches the 
stochastic SIRS model, an active phase with P{I) > 
appears through a continuous transition as b increases. 
In a sufficiently nonlinear regime (large enough a), in- 
creasing b leads to a discontinuous phase transition. The 
nature of the transition can therefore be controlled by a 
single parameter, which is not diffusion. These results 
can be predicted by mean-field analysis and are quali- 
tatively confirmed in simulations of random graphs and 
hypercubic lattices for d > 2. The fact that a discon- 
tinuous transition was not found for d = 1 is consistent 
with previous results for two-state systems with little or 
no diffusion [lS[Il[l3,[2l,[2|. 

We have characterized discontinuous transitions by two 
complementary criteria: first, hysteresis cycles were ob- 
tained and their width estimated by extrapolation for an 
infinite number of Monte Carlo time steps; then, bistabil- 
ity was explicitly confirmed by checking that the system 
trajectory exhibits dependence on the initial conditions. 
In the case of continuous transitions, the width of the 
hysteresis cycles scaled with the ramp rate according to 
recent predictions by Takcuchi p7j . 

Finally, we recall that the exponential coupling in Eq.|4] 
was inspired in the model of Wood et al. Q . While their 
nonlinear phase-coupled stochastic oscillators do undergo 
a phase transition into a synchronized state, we did not 
find sustained collective oscillations with a similar non- 
linearity among pulse-coupled excitable elements. How- 
ever, we did find (albeit in a small parameter region) 
unstable global oscillations and collective excitability in 
the mean-field equations. Simulations of the complete 
graph revealed stochastic oscillations in single runs when- 
ever the active phase corresponded to a stable spiral in 
the mean- field equations. It remains to be investigated 
whether collective excitability and stochastic oscillations 
remain in regular lattices or appear in the transition to 
a small-world regime. 
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